library(ggplot2)

gibbs_betabin <- function(n, a, b, p=0.5, iter=1000){
  x <- matrix(0, iter, 2)
  
  for(k in 1:iter){
    y <- rbinom(1, size=n, p=p)
    p <- rbeta(1, y+a, n+b-y)
    
    x[k,] <- c(y, p)
  }
  
  x
}

set.seed(123)

sp <- data.frame(gibbs_betabin(20, 5, 5))

ggplot(data.frame(Y=sp$X1), aes(Y)) +
  geom_bar(width = 0.5, fill='blue') +
  ylab('Frequency') +
  theme(text=element_text(size=10))
